home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / roots / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  9.6 KB  |  310 lines
  1. /* roots/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_test.h>
  24. #include <gsl/gsl_roots.h>
  25. #include <gsl/gsl_errno.h>
  26. #include <gsl/gsl_ieee_utils.h>
  27.  
  28. #include "roots.h"
  29. #include "test.h"
  30.  
  31. /* stopping parameters */
  32. const double EPSREL = (10 * GSL_DBL_EPSILON);
  33. const double EPSABS = (10 * GSL_DBL_EPSILON);
  34. const unsigned int MAX_ITERATIONS = 150;
  35.  
  36. void my_error_handler (const char *reason, const char *file,
  37.                int line, int err);
  38.  
  39. #define WITHIN_TOL(a, b, epsrel, epsabs) \
  40.  ((fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs)))
  41.  
  42. int
  43. main (void)
  44. {
  45.   gsl_function F_sin, F_cos, F_func1, F_func2, F_func3, F_func4,
  46.     F_func5, F_func6;
  47.   
  48.   gsl_function_fdf FDF_sin, FDF_cos, FDF_func1, FDF_func2, FDF_func3, FDF_func4,
  49.     FDF_func5, FDF_func6;
  50.  
  51.   const gsl_root_fsolver_type * fsolver[4] ;
  52.   const gsl_root_fdfsolver_type * fdfsolver[4] ;
  53.  
  54.   const gsl_root_fsolver_type ** T;
  55.   const gsl_root_fdfsolver_type ** S;
  56.  
  57.   gsl_ieee_env_setup();
  58.  
  59.   fsolver[0] = gsl_root_fsolver_bisection;
  60.   fsolver[1] = gsl_root_fsolver_brent;
  61.   fsolver[2] = gsl_root_fsolver_falsepos;
  62.   fsolver[3] = 0;
  63.  
  64.   fdfsolver[0] = gsl_root_fdfsolver_newton;
  65.   fdfsolver[1] = gsl_root_fdfsolver_secant;
  66.   fdfsolver[2] = gsl_root_fdfsolver_steffenson;
  67.   fdfsolver[3] = 0;
  68.  
  69.   F_sin = create_function (sin_f) ;
  70.   F_cos = create_function (cos_f) ; 
  71.   F_func1 = create_function (func1) ;
  72.   F_func2 = create_function (func2) ;
  73.   F_func3 = create_function (func3) ;
  74.   F_func4 = create_function (func4) ;
  75.   F_func5 = create_function (func5) ;
  76.   F_func6 = create_function (func6) ;
  77.  
  78.   FDF_sin = create_fdf (sin_f, sin_df, sin_fdf) ;
  79.   FDF_cos = create_fdf (cos_f, cos_df, cos_fdf) ;
  80.   FDF_func1 = create_fdf (func1, func1_df, func1_fdf) ;
  81.   FDF_func2 = create_fdf (func2, func2_df, func2_fdf) ;
  82.   FDF_func3 = create_fdf (func3, func3_df, func3_fdf) ;
  83.   FDF_func4 = create_fdf (func4, func4_df, func4_fdf) ;
  84.   FDF_func5 = create_fdf (func5, func5_df, func5_fdf) ;
  85.   FDF_func6 = create_fdf (func6, func6_df, func6_fdf) ;
  86.  
  87.   gsl_set_error_handler (&my_error_handler);
  88.  
  89.   for (T = fsolver ; *T != 0 ; T++)
  90.     {
  91.       test_f (*T, "sin(x) [3, 4]", &F_sin, 3.0, 4.0, M_PI);
  92.       test_f (*T, "sin(x) [-4, -3]", &F_sin, -4.0, -3.0, -M_PI);
  93.       test_f (*T, "sin(x) [-1/3, 1]", &F_sin, -1.0 / 3.0, 1.0, 0.0);
  94.       test_f (*T, "cos(x) [0, 3]", &F_cos, 0.0, 3.0, M_PI / 2.0);
  95.       test_f (*T, "cos(x) [-3, 0]", &F_cos, -3.0, 0.0, -M_PI / 2.0);
  96.       test_f (*T, "x^20 - 1 [0.1, 2]", &F_func1, 0.1, 2.0, 1.0);
  97.       test_f (*T, "sqrt(|x|)*sgn(x)", &F_func2, -1.0 / 3.0, 1.0, 0.0);
  98.       test_f (*T, "x^2 - 1e-8 [0, 1]", &F_func3, 0.0, 1.0, sqrt (1e-8));
  99.       test_f (*T, "x exp(-x) [-1/3, 2]", &F_func4, -1.0 / 3.0, 2.0, 0.0);
  100.       test_f (*T, "(x - 1)^7 [0.9995, 1.0002]", &F_func6, 0.9995, 1.0002, 1.0);
  101.       
  102.       test_f_e (*T, "invalid range check [4, 0]", &F_sin, 4.0, 0.0, M_PI);
  103.       test_f_e (*T, "invalid range check [1, 1]", &F_sin, 1.0, 1.0, M_PI);
  104.       test_f_e (*T, "invalid range check [0.1, 0.2]", &F_sin, 0.1, 0.2, M_PI);
  105.     }
  106.  
  107.   for (S = fdfsolver ; *S != 0 ; S++)
  108.     {
  109.       test_fdf (*S,"sin(x) {3.4}", &FDF_sin, 3.4, M_PI);
  110.       test_fdf (*S,"sin(x) {-3.3}", &FDF_sin, -3.3, -M_PI);
  111.       test_fdf (*S,"sin(x) {0.5}", &FDF_sin, 0.5, 0.0);
  112.       test_fdf (*S,"cos(x) {0.6}", &FDF_cos, 0.6, M_PI / 2.0);
  113.       test_fdf (*S,"cos(x) {-2.5}", &FDF_cos, -2.5, -M_PI / 2.0);
  114.       test_fdf (*S,"x^{20} - 1 {0.9}", &FDF_func1, 0.9, 1.0);
  115.       test_fdf (*S,"x^{20} - 1 {1.1}", &FDF_func1, 1.1, 1.0);
  116.       test_fdf (*S,"sqrt(|x|)*sgn(x) {1.001}", &FDF_func2, 0.001, 0.0);
  117.       test_fdf (*S,"x^2 - 1e-8 {1}", &FDF_func3, 1.0, sqrt (1e-8));
  118.       test_fdf (*S,"x exp(-x) {-2}", &FDF_func4, -2.0, 0.0);
  119.       test_fdf_e (*S,"max iterations x -> +Inf, x exp(-x) {2}", &FDF_func4, 2.0, 0.0);
  120.       test_fdf_e (*S,"max iterations x -> -Inf, 1/(1 + exp(-x)) {0}", &FDF_func5, 0.0, 0.0);
  121.     }
  122.  
  123.   test_fdf (gsl_root_fdfsolver_steffenson,
  124.         "(x - 1)^7 {0.9}", &FDF_func6, 0.9, 1.0);    
  125.  
  126.   /* now summarize the results */
  127.  
  128.   exit (gsl_test_summary ());
  129. }
  130.  
  131.  
  132. /* Using gsl_root_bisection, find the root of the function pointed to by f,
  133.    using the interval [lower_bound, upper_bound]. Check if f succeeded and
  134.    that it was accurate enough. */
  135.  
  136. void
  137. test_f (const gsl_root_fsolver_type * T, const char * description, gsl_function *f,
  138.     double lower_bound, double upper_bound, double correct_root)
  139. {
  140.   int status;
  141.   size_t iterations = 0;
  142.   double r, a, b;
  143.   double x_lower, x_upper;
  144.   gsl_root_fsolver * s;
  145.  
  146.   x_lower = lower_bound;
  147.   x_upper = upper_bound;
  148.  
  149.   s = gsl_root_fsolver_alloc(T);
  150.   gsl_root_fsolver_set(s, f, x_lower, x_upper) ;
  151.   
  152.   do 
  153.     {
  154.       iterations++ ;
  155.  
  156.       gsl_root_fsolver_iterate (s);
  157.  
  158.       r = gsl_root_fsolver_root(s);
  159.  
  160.       a = gsl_root_fsolver_x_lower(s);
  161.       b = gsl_root_fsolver_x_upper(s);
  162.       
  163.       if (a > b)
  164.     gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b);
  165.  
  166.       if (r < a || r > b)
  167.     gsl_test (GSL_FAILURE, "r lies outside interval %g (%g,%g)", r, a, b);
  168.  
  169.       status = gsl_root_test_interval (a,b, EPSABS, EPSREL);
  170.     }
  171.   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
  172.  
  173.  
  174.   gsl_test (status, "%s, %s (%g obs vs %g expected) ", 
  175.         gsl_root_fsolver_name(s), description, 
  176.         gsl_root_fsolver_root(s), correct_root);
  177.  
  178.   if (iterations == MAX_ITERATIONS)
  179.     {
  180.       gsl_test (GSL_FAILURE, "exceeded maximum number of iterations");
  181.     }
  182.  
  183.   /* check the validity of the returned result */
  184.  
  185.   if (!WITHIN_TOL (r, correct_root, EPSREL, EPSABS))
  186.     {
  187.       gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)", 
  188.         r, correct_root);
  189.  
  190.     }
  191. }
  192.  
  193. void
  194. test_f_e (const gsl_root_fsolver_type * T, 
  195.       const char * description, gsl_function *f,
  196.       double lower_bound, double upper_bound, double correct_root)
  197. {
  198.   int status;
  199.   size_t iterations = 0;
  200.   double x_lower, x_upper;
  201.   gsl_root_fsolver * s;
  202.  
  203.   x_lower = lower_bound;
  204.   x_upper = upper_bound;
  205.  
  206.   s = gsl_root_fsolver_alloc(T);
  207.   status = gsl_root_fsolver_set(s, f, x_lower, x_upper) ;
  208.  
  209.   gsl_test (status != GSL_EINVAL, "%s (set), %s", T->name, description);
  210.  
  211.   if (status == GSL_EINVAL) 
  212.     {
  213.       return ;
  214.     }
  215.  
  216.   do 
  217.     {
  218.       iterations++ ;
  219.       gsl_root_fsolver_iterate (s);
  220.       x_lower = gsl_root_fsolver_x_lower(s);
  221.       x_upper = gsl_root_fsolver_x_lower(s);
  222.       status = gsl_root_test_interval (x_lower, x_upper, 
  223.                       EPSABS, EPSREL);
  224.     }
  225.   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
  226.  
  227.   gsl_test (!status, "%s, %s", gsl_root_fsolver_name(s), description, 
  228.         gsl_root_fsolver_root(s) - correct_root);
  229.  
  230. }
  231.  
  232. void
  233. test_fdf (const gsl_root_fdfsolver_type * T, const char * description, 
  234.     gsl_function_fdf *fdf, double root, double correct_root)
  235. {
  236.   int status;
  237.   size_t iterations = 0;
  238.   double prev = 0 ;
  239.  
  240.   gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T);
  241.   gsl_root_fdfsolver_set (s, fdf, root) ;
  242.  
  243.   do 
  244.     {
  245.       iterations++ ;
  246.       prev = gsl_root_fdfsolver_root(s);
  247.       gsl_root_fdfsolver_iterate (s);
  248.       status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev, 
  249.                    EPSABS, EPSREL);
  250.     }
  251.   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
  252.  
  253.   gsl_test (status, "%s, %s (%g obs vs %g expected) ", 
  254.         gsl_root_fdfsolver_name(s), description, 
  255.         gsl_root_fdfsolver_root(s), correct_root);
  256.  
  257.   if (iterations == MAX_ITERATIONS)
  258.     {
  259.       gsl_test (GSL_FAILURE, "exceeded maximum number of iterations");
  260.     }
  261.  
  262.   /* check the validity of the returned result */
  263.  
  264.   if (!WITHIN_TOL (gsl_root_fdfsolver_root(s), correct_root, 
  265.            EPSREL, EPSABS))
  266.     {
  267.       gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)", 
  268.         gsl_root_fdfsolver_root(s), correct_root);
  269.  
  270.     }
  271. }
  272.  
  273. void
  274. test_fdf_e (const gsl_root_fdfsolver_type * T, 
  275.         const char * description, gsl_function_fdf *fdf,
  276.         double root, double correct_root)
  277. {
  278.   int status;
  279.   size_t iterations = 0;
  280.   double prev = 0 ;
  281.  
  282.   gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T);
  283.   status = gsl_root_fdfsolver_set (s, fdf, root) ;
  284.  
  285.   gsl_test (status, "%s (set), %s", T->name, description);
  286.  
  287.   do 
  288.     {
  289.       iterations++ ;
  290.       prev = gsl_root_fdfsolver_root(s);
  291.       gsl_root_fdfsolver_iterate (s);
  292.       status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev, 
  293.                    EPSABS, EPSREL);
  294.     }
  295.   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
  296.  
  297.   gsl_test (!status, "%s, %s", gsl_root_fdfsolver_name(s), 
  298.         description, gsl_root_fdfsolver_root(s) - correct_root);
  299. }
  300.  
  301. void
  302. my_error_handler (const char *reason, const char *file, int line, int err)
  303. {
  304.   if (0)
  305.     printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
  306. }
  307.  
  308.  
  309.  
  310.